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We propose a method for the evaluation of magnetic exchange coupUngs based on noncollinear 
spin-density functional calculations. The method employs the second derivative of the total Kohn- 
Sham energy of a single reference state, in contrast to approximations based on Kohn-Sham total 
energy differences. The advantage of our approach is twofold: It provides a physically motivated 
picture of the transition from a low-spin to a high-spin state, and it utilizes a perturbation scheme for 
the evaluation of magnetic exchange couplings. The latter simplifies the way these parameters are 
predicted using first-principles: It avoids the non-trivial search for different spin-states that needs 
to be carried out in energy difference methods and it opens the possibility of "black-boxifying" 
the extraction of exchange coupUngs from density functional theory calculations. We present proof 
of concept calculations of magnetic exchange couplings in the H-He-H model system and in an 
oxovanadium bimetallic complex where the results can be intuitively rationalized. 



I. INTRODUCTION 

Empirical models based on the Heisenberg spin Hamil- 
tonian are routinely utilized to describe the behavior 
of a variety of magnetic systems. In most cases, these 
simple models are found to fit the experimental data 
very well, provided that the parameters in the model 
Hamiltonian are chosen properly. The set of param- 
eters can include both, external parameters (tempera- 
ture, applied magnetic field, etc.), and internal parame- 
ters (magnetic exchange couplings, magnetic anisotropy, 
etc.). Internal parameters for a particular system can be 
obtained either by fitting experimental data or from first- 
principles electronic structure calculations by mapping 
total electronic energies to the energies of the Heisenberg 
spin Hamiltonian. ^^^^^ In particular, magnetic exchange 
couplings, J, can be obtained considering the isotropic 
Heisenberg Hamiltonian 

^ = -2 ^ J^J S, • S, , (1) 

<ij> 

where and Sj are the (localized) spin operators asso- 
ciated to each magnetic center. 

Perhaps the one of the most interesting manifestation 
of magnetism at the molecular scale can be found in com- 
plexes containing transition metal atoms. Many applica- 
tions have been suggested exploiting these molecular-size 
magnets, such as quantum computation units and high- 
density data storage.^ Due to the relatively large size 
of most complexes of interest, density functional theory 
(DFT)^'^ offers the most efficient alternative for model- 
ing the electronic structure of these systems from first- 
principles. ^^^^^ 

Several approaches had been proposed to extract J 
couplings from DFT energies. According to the spin- 
projected (SP) approach, ^^^^^^ the energies of a two- 
center complex A and B can be related to the J coupling 
as 

Els — Ehs = ^SaSbJab , (2) 



while in the non-projected (NP) approach^^, the energies 
of a two-center complex Sa and Sb can be related to the 
J coupling as 

Els — Ehs = {^SaSb + '^Sb)Jab , (3) 

where Sb < Sa- In Eqs. ([2]) and ©, £; HS is the en- 
ergy of the high-spin state and Els is the energy of 
the low-spin (broken-symmetry) state. Eqs. ([2]) and (|3]) 
can be straightforwardly generalized to a set of equa- 
tions for complexes with multiple magnetic centers. ^-li^ 
While the SP and NP methods are fairly popular, other 
methods have been proposed in the literature such as 
Nishino's approacf>i^, the constrained-DFT approach of 
Rudra et al.^ii^, the Slater's transition state method of 
Dai and WhangboJi and the local spin method of Clark 
and Davidson J^ii^ All these approaches rely on the eval- 
uation of the energy difference between two (for the sim- 
plest case of a bimetallic complex) or more states. The 
evaluation of this energy difference is commonly done by 
carrying out several self-consistent field calculations, one 
for each different magnetic configuration. However, in 
many cases, converging to the right target state could be 
cumbersome, specially for systems containing multiple 
centers with many magnetic configurations. Therefore, 
developing an approach that can be used in a "black- 
box" manner is crucial to systematically explore a large 
set of complexes or complexes containing many magnetic 
centers. 

In this work, we present an approach for the evaluation 
of magnetic exchange couplings based on noncollinear 
spin density functional calculations that allows, in anal- 
ogy to response properties, to express the magnetic ex- 
change couplings as a derivative of the total electronic 
energy of one single state with respect to an external pa- 
rameter, opening the possibility of "black-boxifying" the 
extraction of magnetic exchange couplings from density 
functional theory calculations. 
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II. THEORY AND IMPLEMENTATION 

A. Exchange Couplings as Energy Derivatives 

Let us consider the effective interaction energy between 
two magnetic centers A and B given by the isotropic clas- 
sical Heisenberg model (obtained by taking the expecta- 
tion value of Eq. ([1])), 



= -2Jab SaSbcosO , 



(4) 



where and are the (perfectly localized) magnetic 
moment vectors, Jab is the exchange coupling constant, 
and is the angle between and Sb- From Eq. (|4]), 
one can trivially obtain Jab from the second derivative 
of Eab with respect to at the equilibrium points. 



J AE 



2SaSb 
1 



2SaSb V dO^ 



(5) 



6>=180° 



These simple relations provide a direct path to the eval- 
uation of magnetic exchange couplings Jab from density 
functional calculations if Eab in Eq. ([5]) is replaced by 
the total Kohn-Sham (KS)^^ energy of the system, E^^. 
Therefore, assuming that the electronic system depends 
on as an ideal Heisenberg model (the validity of this 
assumption will be discussed in the next Section), one 
can express the exchange coupling constant Jab in terms 
of an energy derivative as 



Ja 



or 



Ja 



2SaSb V dO'^ 



1 fd^E""' 
'2SaSb V dO^ 



(6) 



(7) 



6'=180° 



where the angle 6 in the DFT framework is defined as the 
angle between the local magnetization vectors and 
S^. Another related method based on the Green's func- 
tion formalism for crystals has been proposed by Liecht- 
enstein et al^ 



The two-component spinors introduce the freedom in the 
spin-dependence of the KS system that allows for local 
rotations of the spin density characterized by 6> 7^ and 
6 ^ 180°, i.e. noncollinear spin densities.— i^i^i^i^ The 
local magnetization vectors and S b can be written as 



where 



Sa,b = j 6^VWA,B(r)s(r), 



s(r) = ^\{v)cT^,(v) 



(10) 



(11) 



iGocc 



is the spin-density vector and Wa,b{'^) is a scalar weight 
function that determines each local magnetic site. It is 
important to recall that the magnetic centers A and B 
represent a group of one or more atoms. In Eq. (pTj) . 
a = ((Ja^, (jy^cTz) is the 2x2 Pauli matrices vector. 

Having defined the local magnetic moments, the sec- 
ond step is to find the dependence of the total electronic 
energy upon local rotations of the spin density. This is 
done by constraining the direction of the local magneti- 
zations Sa and by means of the Lagrange multipli- 
ers technique. To this end, we construct the Lagrangian 
functional A 



A[{*(r)},AA,AB]=i?-^[{f(r)}]- 
■ (S^ X z) - As ■ (Sb x e^) , 



(12) 



where e^* = sin ^ x + cos 6> z is a unity vector to which 

is constraint to be parallel to, and for simplicity 
has been chosen to be constraint to the z direction, as 
schematized in Fig. [TJ Here < ^ < 180° is consid- 
ered as an external parameter for which {dE/dO)o^o = 
{dE/d6)e=i80° = 0. Eq. ([12]) can be readily generalized 
for the case of many magnetic centers and arbitrary unity 
vector directions. For the case of two magnetic centers 
(and for the purpose of this work), Eq. (p!2]) does not im- 
ply any loss of generality. In Eq. ([12]), £^^^[{^(r)}] rep- 
resents the KS energy of the system which is, in practice, 
a functional of the set of occupied KS orbitals {^(r)}. 
Stationary solutions of A for a given (fixed) imply 



dA 
dXA 



= Sa X z = . 



(13) 



B. Constraint Noncollinear Spin-DFT Calculations 

To evaluate the dependence of E^^ on 6>, we first in- 
troduce two-component spinors as Kohn-Sham orbitals. 



*i(r) 



^J(r) 
V'l(r) 



(8) 



where tpj (r) and ?/;| (r) are spatial orbitals expanded in a 
linear combination of atomic orbitals, 



V'r(r) = ^c-(/.^(r) {u; =14). 



(9) 



and 



^A 



dA 
dXn 



Sb X = , 



SE'' 



(14) 



WA(r)AA • (o- X z) 



Wb{v)Xb • {(T X e^)] ^^(r) = (i G occ). (15) 

While Eqs. ([T3|) and restore the constraint condi- 
tions, Eq. ([T5]) combined with the orthonormality condi- 
tion for the spinors yields a modified set of KS equations 



3 



(in terms of two-component spinors) that include the two 
additional terms inside the square brackets on the left- 
hand side of Eq. (p!5|) . 



T + Fat + J + Kcc - WA{r)XA ■ {a x z) - 

WbWAb • {a X e^)l ^,(r) = e,^,(r) , (16) 



where T = — 1/2V^ is the kinetic energy, Vn is the 
electron-nuclei potential, J is the Coulomb (or Hartree) 
potential, and Vxc is the exchange-correlation (XC) po- 
tential. The sum of the first four terms inside the square 
brackets is the standard KS Hamiltonian, while the two 
additional terms can be interpreted as a potential orig- 
inated in a torque exerted on the local magnetic mo- 
ments and S^. It should be noted that other ap- 
proaches had been proposed in the literature to con- 
straint the direction of the local mangetization in spin 
DFT calculations.27^28,29 

It is important to note that since the constraint condi- 
tions are linear in the spin density vectors, the additional 
terms in Eq. (p!6|) depend implicitly on the orbitals only 
through Xa and A^, which simplifies the implementa- 
tion. 




FIG. 1: Schematic representation of the constraint vectors 
employed for the local rotations of the spin density. 



Neglecting spin-orbit interaction, T, Vat, and J are 
diagonal in the 2x2 spin space, and thus the only 
term in Eq. (p!6|) that couples ipj and ip^ is Vxc In a 
previous work, we have generalized Vxc for noncollinear 
magnetizations^ assuming that the XC energy depends 
on the local variables in the same manner as in the 
standard collinear (spin- unrestricted) case, and impos- 
ing the condition for the XC energy to be invariant un- 
der rigid rotations of the spin density. In that work, we 
have derived Vxc for general energy functionals contain- 
ing a variety of ingredients beyond the local-spin den- 
sity approximation (LSDA) and the generalized-gradient 
approximation (GGA), such as meta-GGAs and hy- 
brid density functionals. This same generalization is 
adopted throughout this work. Other implementations 
based on plane- wave a^^i^^i^^i^^ and Gaussian-type or- 
bitals can be found in the literature for a variety of 
applications ^^i^i^i^i^ii^ 

The constraint vectors can be chosen without loss of 
generality to lay in the x — z plane. Hence, as spin-orbit 



interaction is not included in the Hamiltonian, the two- 
component spinors are purely real. As a consequence, in 
this scheme the orbital magnetization of the solutions is 
always zero. 

The value of A at the stationary solutions for a fixed 
given by Eqs. (p!3|) - (p!5|) can be directly associated to 
the energy of the KS system in the presence of the con- 
straints, E^^{6). The only arbitrariness in the formula- 
tion is the choice of the weight function Wa,b{^) and the 
fragments A and B. Since in practice it is necessary to 
evaluate the matrix elements of Wa,b{^) in the atomic 
orbitals basis set ^^(r), it is convenient to define and 
Sb using population analysis. 



(17) 



where P^^y is the spin-density matrix vector whose Carte- 
sian components are 



and 



px _ pTi I pit 
py _ j(pU _ piTA 



pz _ pTT I pii 



(18) 
(19) 

(20) 



written in terms of the 2x2 density matrix 

^;"'= E {Lo,j=\,i). (21) 



iGOCC 



For this work, we employ Lowdin population analysis*^ 
From the expression for the atomic spin-population given 
by Eq. (pT|). {y\^A,B)^v can be obtained as 



(Wa,b)mi^ 



{t] = x,y OT z) . (22) 



Thus, using the Lowdin partitioning, it is straightfor- 
ward to show that the matrix elements of Wa,b{'^) can 
be calculated from the atomic orbitals overlap matrix 
{S)xa = J d^rM^)M^) as 

{Wa,b),u= ^ (§^/^)am(§'/').a. (23) 
xeA,B 

It is worth commenting that atomic spin-densities are 
less sensitive to the choice of the population method than 
atomic densities. Since only the direction of the atomic 
spin-density is relevant in our method for the calcula- 
tion of magnetic exchange couplings, we expect even less 
sensitivity with the coice of the population method. 

The set of spinors that simultaneously satisfy Eq. ([16]) 
and Eqs. (p!3|) and (p!4|) needs to be determined self- 
consistently since J and Vxc depend on the spinors. To 
obtain the KS Hamiltonian, we add the additional con- 
straint terms of Eq. ([16]) to the standard KS Hamilto- 
nian, initially using a guess for the Lagrange multipliers 
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and A b • Then we determine the optimal A^ and A b 
such that Eqs. (p!3|) and ([E]) are satisfied, using the den- 
sity matrix obtained from diagonahzing the KS Hamil- 
tonian of the constraint system to evaluate and Sb- 
This is carried out using the steepest descent method to 
find the minimum of A as a function of Xa and Xb- Once 
the optimal values of the Lagrange multipliers are deter- 
mined, and provided that self-consistency is not achieved, 
we proceed to the next self-consistent iteration using the 
density matrix from the previous iteration. The process 
stops once the criteria for changes in the density matrix 
and total energy are met. Several consistency checks were 
performed to verify the robustness of our code. We have 
implemented this scheme in the Gaussian Development 
Version program. 

Using this methodology, E^^iO) was calculated for 
small values of 9 around 6> = and = 180° and magnetic 
exchange couplings were obtained from the quadratic co- 
efficient of a polynomial fit. It is important at this point 
to mention that there are cases where existing approxi- 
mate density functional methods have difficulties in rep- 
resenting the LS state (0 = 180°)^^^ The Kohn-Sham 
determinant in these cases correspond to a "broken- 
symmetry" solution that mixes two or more eigenfunc- 
tions of the S'^ operator. However, for the HS state 
{6 = 0), it is customary accepted that approximate den- 
sity functionals provide a reliable representation. Thus, 
even though for comparison purposes in the next Section 
we show our results using both the HS and LS states 
as reference, for the practical extraction of magnetic ex- 
change couplings in the DFT framework, this method is 
expected to work more reliably using only the HS as the 
reference state. 



III. PROOF OF CONCEPT CALCULATIONS 

We first tested our methodology in the H-He-H linear 
model system with a distance H-He of 1.625 A, consid- 
ering the outer H atoms as magnetic centers A and B 
{Sa = 5'b = 1/2). In Table [J we show our results for the 
magnetic exchange couplings calculated from d?E^^ / dO'^ 
using the 6> = {{sj) = 1) state and the 6 = 180° 
{{Sz) = 0) state, J^^ and J^^, respectively. All cal- 
culations were carried out with the 6-31 IG** Gaussian 
basis set.^^ For comparison, in Table [J we include results 
for the LSDA (Dirac exchange plus the parametrization 
of Wosko, Wilk, and Nusair^^ for correlation), the BLYP 
realization of the GGA (Becke's 1988 functional^^ for ex- 
change and the correlation functional of Lee, Yang, and 
Parri^), and for the mhYV^^^^ hybrid functional. 
For all functionals, exchange couplings calculated from 
the energy derivatives, J^^ and J^^, are in very close 
agreement to the exchange coupling calculated from the 
energy difference, J^^. The difference can be attributed 
to both, the intrinsic accuracy of the numerical differ- 
entiation method and to the fact that J^^ and J^^ are 
expected to be identical to J^^ only in the case where 



DFT describes the electronic system as an ideal Heisen- 
berg model. The small discrepancy between J^^, J^^ 
and J^^ can be understood in terms of the localized na- 
ture of the magnetization on the H atoms in this model 
system and provides a measure of how well the electronic 
system mimics the behavior of an ideal Heisenberg model. 

In Table [J we report J^^ as calculated from the SP 
formula, Eq. ([2]), since it offers a direct comparison be- 
tween J^^ and J^^, and J^^. It is worth mentioning 
that in the ideal case of a perfect Heisenberg system, 
A^; = E{e = 180°) - E{e = O) is related to J^^ and J^^ 
according to 

Therefore, dPE/dO'^ provides a measure of /S.E that 
can be evaluated without explicitly converging the self- 
consistent procedure to the LS state. 



TABLE I: Magnetic exchange couplings (in meV) calculated 
from energy derivatives and from energy differences for the 
H-He-H system. 







LSDA 


BLYP 


B3LYP 


JHS ^ 


1 (d^E\ 
2SaSb \~dO^Je^O 


-95.8 


-74.0 


-60.8 


JLS ^ 


1 (d''E\ 
2SaSb \~dO^ J 0^180° 


-101.7 


-76.6 


-61.2 


jAE ^ 


4SaSb 


-99.8 


-76.9 


-63.5 



Our second proof of concept was carried out in the ox- 
ovanadium(IV) dimer [(/i-OCH3)VO(ma)]2. This com- 
plex shows a strong antiferromagnetic coupling of about 
— 13.3 meV, as measured by temperature-dependent 
magnetic susceptibility experiments.^^ Here we employed 
Ahlrich's triple-zeta valence basis set for for the V 
atoms and Ahlrich's double-zeta valence basis for first- 
row atom a^^i^^ , as obtained from Ref. 48. This basis 
was shown to provide reliable results in practical cal- 
culations of exchange couplings.— Atomic coordinates 
were taken from experimental crystallographic data4^ In 
Figs. [2] and [3] we present a plot of E^^ as a function 
of {0 < < 180°) for LSDA and B3LYP, respec- 
tively. In both figures, E^^{6) follows closely a cosine 
function connecting the HS and LS extrema, indicating 
that both functionals capture the Heisenberg behavior 
of the oxovanadium complex. Related investigations in 
periodic systems using the LSDA can be found in the 
literature^^i^ 

For the plots shown in Fig. [2] and Fig. [3] we have chosen 
as magnetic centers A and B {Sa = Sb = 1/2) both sets 
of V and apical O atoms since most of the spin density 
in this complex is localized on this moiety, as shown in 
Fig.m However, it is worth to mention that by choosing 
the metal atoms only as magnetic centers A and B the 
changes in the plots are unappreciable and the calculated 
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9 (degrees) 



TABLE II: Magnetic exchange couplings (in meV) calcu- 
lated from energy derivatives and from energy differences for 
the vanadium bimetallic complex. The corresponding exper- 
imental value is —13.3 meV." 





LSDA 


B3LYP 


jHS 1 / d^E \ 

2SaSb \~dO^Je^O 
jLS _ if d'^E\ 

2SaSb V"^/ 0=180° 
jAE E{e^lSQ°)-E{e^Q) 


-46.6 
-41.9 
-44.1 


-11.4 
-11.2 
-11.3 


^Taken from Ref . ^5. 



only 2.75° (for = 90°). 



FIG. 2: LSDA energy change as a function of the angle be- 
tween the local magnetic moments obtained from a constraint 
noncollinear spin density functional calculation in the oxo- 
vanadium complex (Fig. |4|. The solid line shows the (ideal) 
cosine function connecting the AF and FM extrema. 




-12 ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

20 40 60 80 100 120 140 160 180 

9 (degrees) 



FIG. 3: Same as Fig.[2]for B3LYP. 



magnetic exchange couplings J^^ and J^^ vary very lit- 
tle. For instance, for LSDA, magnetic exchange couplings 
change (in meV) from J^^ = —46.6 to J^^ = —46.9 and 
from J^^ -41.9 to J^^ -42.2 when using the V 
and O atoms or the V atoms only, respectively. This 
indicates that our method is not sensitive upon a par- 
ticular choice of the magnetic centers and shows, in this 
sense, robustness. One physical explanation for this fact 
is that the spin polarization of the light atoms surround- 
ing a metal atom is "dragged" by the strong magnetic 
coupling with the neighbor metal center and therefore, it 
tends to align parallel (or antiparallel) to the magneti- 
zation of the metal center. For the case in study, if the 
constraint is applied on the V atoms only, the angle of 
the spin polarization on the apical O atom deviates from 
the direction of the constrain vectors by a maximum of 




FIG. 4: Spin density isosurface of the HS (a) and LS (b) 
states of the oxovanadium complex (Fig.|4|. Red corresponds 
to t and blue to |. The isosurface represents a spin density 
of 0.01 (Bohr~^). Magnetic centers A and B chosen for the 
calculations are also indicated. 

A careful comparison of Fig. [2] and Fig. [3] evidences a 
larger deviation from the ideal cosine function for LSDA 
than for B3LYP. The largest differences from the cosine 
function are approximately 0.6 meV and 0.04 meV for 
LSDA and B3LYP, respectively, and occurs for = 90° 
in both cases. This is not surprising since LSDA yields 
electron (total and spin) densities more delocalized than 
B3LYP (the Lowdin atomic magnetic moments at the V 
atoms for the LS state are 1.00 /i^ and 1.10 /i^ for LSDA 
and B3LYP, respectively) and therefore, one can expect 
that the B3LYP energy follows the Heisenberg behavior 
more closely than its LSDA counterpart. Localization of 
the spin-density also reduces the calculated magnetic ex- 
change couplings, as shown by Martin^^, Ruiz^^^^^, and 
demonstrated by Rudra et at. by explicitly constraining 
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the local magnetization of the LS state. As shown in 
Table [ni The difference between J^^ and J^^ is 4.7 meV 
for LSDA, while it is only 0.2 meV for B3LYP. Thus, in 
contrast to the perfectly localized case, a more delocal- 
ized magnetization yields to larger deviations from the 
ideal Heisenberg behavior and hence greater differences 
between J^^ and J^^ and at the same time to larger 
exchange couplings. 

IV. CONCLUSIONS 

We have proposed a method for the calculation of mag- 
netic exchange couplings from noncollinear spin density 
functional calculations that employs the second deriva- 
tive of the electronic energy of a single state with respect 
to a parameter, Eqs. (|6]) and ([7j). Within this approach 
there is no need to search for different self-consistent so- 
lutions of spin-states as it is commonly done in methods 
based on energy differences, such as the SP or NP meth- 
ods, Eqs. ([2]) and ([3|). Our method utilizes perturbation 
theory for the evaluation of magnetic exchange couplings 
and therefore, in combination with standard analytic sec- 
ond derivatives techniques, it can potentially be used to 



compute exchange couplings very efficiently, opening the 
possibility of "black-boxifying" the extraction of mag- 
netic exchange couplings from density functional theory 
calculations. 

Our proof of concept calculations show very promis- 
ing results. For the cases studied we found that our 
method reproduces exchange couplings obtained from 
the spin-projected approach based on energy differences. 
As expected from physical grounds, the agreement be- 
tween both methods improves when the DFT description 
of the interaction between the magnetic centers is more 
Heisenberg-like. In this sense, the curve E^^{0) provides 
a quantitative measure of how well the electronic system 
mimics the behavior of an ideal Heisenberg model. 
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